#multilevel analysis
nullmodel <- lmer(happy ~ (1|country),data=wv6,REML=FALSE)
summary(nullmodel)
##empty random intercept
rho(nullmodel)
countries <- unique(wv6$country)
regs <- sapply (
  X = 1:length(countries),
  FUN <- function(i) {
    lm <- lm(happy ~ fam_trust+par_trust+gen_trust+married+health+income+married+male,data=wv6,country==countries[i])
    c(coef(lm)[1:2],mean(lm$fitted))
  })
rownames(regs)=c("b_const","b_fam_trust","mean.fitted")
overall<-round (rowMeans(regs,na.rm=T),2)
overall
wv6$country
## plot the above
regs.df <- data.frame(t(regs))
no_legend<-theme(legend.position="none")
g_fam_trust <- ggplot(regs.df,aes(x=b_fam_trust,y=mean.fitted,color=factor(mean.fitted)))
g_fam_trust<-g_fam_trust+no_legend
g_fam_trust+geom_point()
mean_point<-annotate("point",x=overall[2],y=overall[3],color="navyblue",size=8)
g_fam_trust+geom_point()+mean_point
##plot for par_trust
regs <- sapply (
  X = 1:length(countries),
  FUN <- function(i) {
    lm <- lm(happy ~ par_trust+fam_trust+gen_trust+married+health+income+married+male,data=wv6,country==countries[i])
    c(coef(lm)[1:2],mean(lm$fitted))
  })
rownames(regs)=c("b_const","b_par_trust","mean.fitted")
overall<-round (rowMeans(regs,na.rm=T),2)
overall
regs.df <- data.frame(t(regs))
no_legend<-theme(legend.position="none")
g_par_trust <- ggplot(regs.df,aes(x=b_par_trust,y=mean.fitted,color=factor(mean.fitted)))
g_par_trust<-g_par_trust+no_legend
g_par_trust+geom_point()
mean_point<-annotate("point",x=overall[2],y=overall[3],color="navyblue",size=8)
g_par_trust+geom_point()+mean_point
##plot for gen_trust
regs <- sapply (
  X = 1:length(countries),
  FUN <- function(i) {
    lm <- lm(happy ~ gen_trust+par_trust+fam_trust+gen_trust+married+health,data=wv6,country==countries[i])
    c(coef(lm)[1:2],mean(lm$fitted))
  })
rownames(regs)=c("b_const","b_gen_trust","mean.fitted")
overall<-round (rowMeans(regs,na.rm=T),2)
overall
regs.df <- data.frame(t(regs))
no_legend<-theme(legend.position="none")
g_gen_trust <- ggplot(regs.df,aes(x=b_gen_trust,y=mean.fitted,color=factor(mean.fitted)))
g_gen_trust<-g_gen_trust+no_legend
g_gen_trust+geom_point()
mean_point<-annotate("point",x=overall[2],y=overall[3],color="navyblue",size=8)
g_gen_trust+geom_point()+mean_point
## random slopes
lmer1<- lmer(happy ~ fam_trust+par_trust+gen_trust+married+male+health+income+(1+gen_trust|country),data=wv6,REML=FALSE)
summary(lmer1)
###Graph: correlation between coef of trust and intercept for each country
fam_trust_vs_const <- ggplot (regs.df,aes(x=b_fam_trust,y=b_const))+no_legend
fam_trust_vs_const+geom_point(aes(color=factor(b_const)))+stat_smooth(method="lm",se=F)
par_trust_vs_const <- ggplot (regs.df,aes(x=b_par_trust,y=b_const))+no_legend
par_trust_vs_const+geom_point(aes(color=factor(b_const)))+stat_smooth(method="lm",se=F)
gen_trust_vs_const <- ggplot (regs.df,aes(x=b_gen_trust,y=b_const))+no_legend
gen_trust_vs_const+geom_point(aes(color=factor(b_const)))+stat_smooth(method="lm",se=F)
###Compare random intercept vs. random intercept+slope
anova(lmer1,lmer2)
###Add cross-national interactions
lmer2<- lmer(satlife ~ ln.GDP*par_trust+ln.GDP*fam_trust+married+male+health+income+education+(1+fam_trust+par_trust|country),data=wv6,REML=FALSE)
summary(lmer2)
lm2 <- lm(satlife ~ Gini*agg_trust+married+male+education+income+health,data=wv6)
summary(lm2)

